Statistics of trajectories in two-state master equations 
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We derive a simple expression for the probability of trajectories of a master equation. The 
expression is particularly useful when the number of states is small and permits the calculation of 
observables that can be defined as functionals of whole trajectories. We illustrate the method with a 
two-state master equation, for which we calculate the distribution of the time spent in one state and 
the distribution of the number of transitions, each in a given time interval. These two expressions 
are obtained analytically in terms of modified Bessel functions. 

PACS numbers: 02.50.Ga, 05.10.Gg 



The evolution of many systems in physics, chemistry 
and biology is properly described by master equations. 
This description is adequate when the system under con- 
sideration has discrete states and when the rate of jump- 
ing from one state to another does not depend on the 
history of the system, i.e. when the Markov property 
holds. In recent years, this description has been suc- 
cessfully applied to a plethora of new problems in sev- 
eral fields. As examples, master equations are commonly 
used in biochemistry to understand the fluctuations of 
chemical concentrations inside the cell In statistical 
physics, they can provide a simple description of non- 
equilibrium systems, useful for testing the validity of fluc- 
tuation relations [3, S 3 • 

From a technical point of view, master equations now 
constitute a well established field of research, and many 
techniques have been developed which permit their ana- 
lytical or numerical treatment 0, S 0] ■ I n complicated 
cases, these techniques permit calculation of the steady- 
state probabilities, P n , of being in state n. In simpler 
cases, it is sometimes possible to solve equations in time 
in order to determine the propagator, p(n,i|no,0), that 
gives the probability of being in state n at time T starting 
from a state no at time t = 0. 

For many practical purposes, determination of the 
propagator is sufficient, since many interesting observ- 
ables can be expressed as a function of the propagator. 
There are, however, observables that cannot be obtained 
conveniently from the propagator, including in particular 
quantities which are more easily expressed as functionals 
of entire trajectories. Examples include the distribution 
of the time spent in a given state and the probability 
of observing a given number of transitions, both for a 
fixed time interval. Functional methods are well known 
for continuous stochastic process, where techniques have 
been developed in parallel to those used in quantum me- 
chanics [H . There are fewer examples of functional meth- 
ods for discrete processes [H [l(| . These methods are of- 
ten field theoretic in nature and involve complications 
such as rcnormalization which one would like to avoid in 
simple discrete systems. 

In this paper, we present a simple way to calculate 
probabilities of the trajectories of master equations. The 



method is straightforward, rigorous and does not require 
any specific assumptions on the equation. It is particu- 
larly useful when the number of states available to the 
system is small, where it is possible to obtain closed ana- 
lytical expression for several interesting observables. We 
study as example of our method general two-state sys- 
tems that, despite their simplicity, have many non-trivial 
applications in problems related to single-molecule spec- 
troscop y ( see, e.g. [II], fl2| |) and biophysics (see, e.g., 
[HI, ELlHI)- Specifically, we calculate the probability of 
observing a given number of transitions, N, in a time, T, 
and the distribution of time spent in one of the two states 
in a time T. Each of these quantities can be expressed 
in terms of modified Bessel functions. 
We consider a master equation: 

j t P n = WmnPm ~ W nm P n , (1) 

m 

where P n {t) is the time-dependent probability of being 
state n and W mn is the transition rate from state m to 
n. For convenience we also define: 

WZ* = Y i W nii (2) 

i 

the total out-rate of state n. The probability that, in a 
time T, the trajectory visits a pre-determined sequence 
of states no, n±, ... then becomes 

V(no, m, ... njv; T) = I dt\ I dt 2 ■ ■ ■ An x 

xe-^'w„ „ ie -^( t2 - tl ) . . . WnN _ inN e- w ^ T -^.(i) 
This expression can be understood by noticing that the 
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FIG. 1: Trajectory of a master equation as a function of time. 
The integral in eq. ((3j is over the times of the N transition 
points. 

integrand represents the probability density in time of 
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the N consecutive transitions according to the master 
equation (see Fig. 1). By summing over all trajectories 
having prc-dctcrmincd properties, one can reconstruct 
the full statistics of the stochastic process. An obvious 
example is the propagator, that can be evaluated as the 
sum over all trajectories that start in a given state, no, 
at time t = and end in a state n/ at time T: 

oo 

p(nf,T\no,0) = ^2 ^ V(n ,rii, . . .n N = n/;T) . 

JV=0{n 1 ...n JV -l} 

(4) 

Note that all probabilities are properly normalized; in 
particular, the propagator satisfies the closure relation 
E„,P(»/,r|no,0) = l. 

The above expressions become particularly useful when 
the number of distinct states visited by the system is 
small. In this case, it is convenient to rearrange the inte- 
grals in eq. ([3]) by grouping together all time intervals in 
which the system is in the same state. If ki is the number 
of times the system visits state i on a given trajectory, 
one finds 



V{n ,ni,n 2 . . - n N ;T) = W noni W ni „ 2 . . 
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where the index i runs over all states visited by the sys- 
tem at least once in the given sequence. 

As an example, we consider the simple case of a mas- 
ter equation with two states, + and — , with transition 
rates k + (from — to +) and fc_ (from + to — ). In spite 
of its simplicity, this case is of interest for many physi- 
cal and biological problems [ill, El El HI El ■ We will 
show that eq. (O allows analytic calculation of the prob- 
abilities of different classes of trajectories. This makes 
it possible, for example, to obtain closed expressions for 
the probability of observing a given number of transitions 
in a time T and for the probability of spending a given 
time in states ± during a time T. For convenience, we 
introduce here the total rate kr = /c+ + fc_ and the equi- 
librium probabilities, = k + /kr and Pt q = k_/kx- 
In this case, we can immediately write the probabilities 
of all possible trajectories according to eq. (0). The sim- 
plest trajectories are evidently those in which there is no 
transition in the interval 

[0,T]: 



V(+;T) = c 
V(-,T) = c 



-k-T 



-k+T 



(0) 



slightly different expressions are obtained for N even and 
for TV odd. The result is: 



V(-,T,t+,N even ) 
V(-,T, t+,N odd ) 

V{+,T,t+,N even ) 
V(+,T, t +1 N odd ) 



[fc + (T-t + )]f(fc_t + )(f- 1 ) 
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fc_ e" r . (7) 



where r = [k-t+ + k+(T — t+)]. These equations describe 
all trajectories with N > while cqn. ^ describes the 
two trajectories with N = 0. Note, however, that eqn. (J6]) 
describes probabilities while cqns. arc probability den- 
sities in t+. To obtain consistent notation, the two ex- 
pressions in eq. © should be multiplied by 5(t+ — T) and 
S(t + ), respectively. This formalism allows us to calculate 
the distribution of time spent in a state during a time 
interval T, g(t±\T). Drawing the initial state from the 
equilibrium distribution (Pf 3 , Pl q ), we find 



!(t + \T) = P? £ P(+,t + ,N)+Pl« ^Tv(-,T,t + ,N). 



N=0 



N=0 



(8) 

Inserting eqns. © and ((7|) into this expression and sum- 
ming the series, we obtain 



g(t+\T) = Pl q e- k + T 5(t+ 



t+ 1 (T-i+) J k T 



P; q c- k + T S{T-t+) + 
* Il{2z ) + 2 J?±^I (2z)\ , (9) 



where r is same as in eq. ([?]) ■ 



+ k + (T 



and Iq(z) and Ii(x) are modified Bessel functions. Notice 
that this result can be obtained in a less direct way by 
means of the Anderson formalism [l6[. Note also that 
the propagators can be obtained by an integration over 
t + of the various terms contributing to g{t±\T). 

In Fig. ^ we plot the function g(t+\T) for several val- 
ues of the parameters, and we compare it with simu- 
lations of the master equations. In all cases studied, 
there is perfect agreement between the simulations and 
the present analytic result. 

An interesting limit of eq. @ is that of large T. 
Using the asymptotic expression lim. ! ;_ +00 [/^(a;)] = 
exp [xj , we see that the leading term in 1/T is 



The determination of general trajectories is simplified by 
having only two states, since trajectories can only al- 
ternate between them. It is then convenient to classify 
trajectories according to: a) the initial state (+ or — ), b) 
the total time T, c) the total time spent in state ±, t±, 
and d) the total number of transitions, N. This is suffi- 
cient to characterize a general term in eq. ([S]). Note that 



ff (*+|T) « /A^e-M^VM^) 2 (10) 

V 7TZ kT 

which has exponential tails expected from large deviation 
arguments [17| . 

Another issue that can be addressed in this framework 
is the probability, h(N), of observing precisely N transi- 
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FIG. 2: A comparison of a simulation and eq. (JUJ for the func- 
tion g(t+\T). The parameters are T = 5 (top figures) and 
T = 50 (bottom figures). The rates are fc+ = fc_ = 0.5 (left 
figures) and k- = 0.8 and k+ = 0.2 (right figures). Lines (red 
on-line) are the analytic curves; the black point are averages 
over 10 7 simulations of the master equation. Notice the effect 
of the Delta functions (first and last point) in the top figures, 
where the probabilities are in log scale. (We do not plot the 
delta functions in the analytic curves). 
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FIG. 3: The probabilities, h(N), of observing N transitions 
in a time T = 50. Transition rates are (left) k+ = fc_ = 
0.5 (right) and k+ = 0.2 and fc_ = 0.8 (left). The points 
represent statistics collected over 10 7 simulations. The solid 
lines (red on-line) are the analytic results of eqns. (|I2|I and 
(|13p . The left figure corresponds to the symmetric limit with 
k + + = k- = 0.5 in which both distributions collapse into a 
Poisson distribution. Notice the even-odd asymmetry in the 
right figure. 



tions in a time interval of T. This is simply: 

h(N) = [ dt + [Pl q V(-,T,t + ,N) + Pl q T(+,T,t + ,N)}. 
Jo 

(11) 

Using the above expressions for the various terms and 
performing the integral, we find two different expressions, 
one for N odd: 



((TV- l)/2)!fc T \k--k +i 

xe-^ +k -W 2 I N/2 (0 (12) 

and one for N even: 

h(N) _ v^T(fc_fc+)W 2 ) / T \ L ^ I B -(fc + +fc.)g 
'H J V — 2k T (N/2)\ \k--k + J e 

X [(*_ + ft+) V-l)/2 (0 + (ft- - ft+) I(N+l)/2 (C)] (13) 

with C = (ft- - k + )T/2. 

In evaluating the two expressions above, care must be 
taken to pick up the proper branch of the half-integer 
powers according to the requirement that the function 
h(N) should be real and positive. In Fig. ([3]) we com- 
pare the distribution h(N) with simulations of the mas- 
ter equation. Here, too, perfect agreement is found. The 
left panel shows a symmetric case with k + = k— =0.5 for 
which eqns. (|12p and (fl3|) each have as limit a Poisson dis- 
tribution, h(N) = X N e~X/N\ with A = T/k + = T/k— 
The right panel, for the case k + = 0.2 and fc_ = 0.8, is 
less trivial. The asymmetry in the rates is reflected in 
a difference between the distributions for N even and N 
odd. This corresponds to the physical fact that one of 



the states is short-lived and the other long-lived, so that 
is more likely to observe an even number of transitions. 
In the asymmetric case, accurate numerical studies indi- 
cate that the average number of transitions is N = T/k 
with k = k + k_/[2(k + + k-)]. 

In summary, we have shown that the probability distri- 
butions associated with the trajectories of master equa- 
tions can be expressed in general as a product over single- 
state properties. This can be particularly useful for sys- 
tems composed of a few states as demonstrated by the ex- 
act determination of several statistical quantities of two- 
state master equations for which results can be expressed 
simply in terms of modified Bessel functions. While the 
methods presented here can be applied to the evalua- 
tion of individual trajectories in more complex problems, 
summation over all trajectories becomes increasingly dif- 
ficult as the number of states increases. If, however, al- 
most all rates are small, the dynamics of the system can 
be dominated by a relatively small number of trajecto- 
ries. For example, this is often the case in chemical ki- 
netics, where average reaction paths may be well defined 
even for high-dimensional dynamics [la ] . In such cases, 
our methods could provide a way to detect these domi- 
nant trajectories and to assess their probabilities. 



Acknowledgments 

We would like to thank E. Barkai for pointing out rel- 
evant references. S. P. wishes to thank J. Fcrkinghoff- 
Borg, J. Fonslet, M. H. Jensen and S. Krishna for stim- 
ulating discussion. 



4 



[1] M.B. Elowitz, A.J. Levine, E.D. Siggia, P.S. Swain, Sci- 
ence 297(5584) pp.1183-1186 (2002). 

[2] J. Lebowitz and H. Spohn, Jour. Stat. Phys. 95, pp 333- 
364 (1999). 

[3] M. Esposito, U. Harbola and S. Mukamel, Phys. Rev. E 
76, 031132 (2007). 

[4] J. Harris and G. M. Schutz, J. Stat Mech. P07020 (2007). 

[5] C. W. Gardiner, Handbook of stochastic methods 
(Springer, Berlin 1983). 

[6] D. T. Gillespie, J. Phys. Chem. 81, pp. 2340-2361 (1977). 

[7] S. Pigolotti and A. Vulpiani, Jour. Chem. Phys. 
128(154114), (2008). 

[8] See, e.g.,C. Itzykson, J. M. Drouffe, Statistical Field The- 
ory vol. 1, Cambridge University Press (1989). 

[9] L. Peliti, J. Phys. (Paris) 46, 1469 (1985). 
[10] J. Cardy, in The Mathematical Beauty of Physics, edited 
by J. Drouffe and J. B. Zuber (World Scientific, Singa- 



pore 1997) 

[11] G. Margolin and E. Barkai, Phys. Rev. Lett. 94, 080601 
(2005). 

[12] F. Shikerman and E. Barkai, Phys. Rev. Lett. 99. 208302 
(2007). 

[13] G. Bonnet, O. Krichevsky, and A. Libchaber, Proc. Natl. 

Acad. Sci. 95(15) pp. 8602-8606 (1998). 
[14] R. Zwanzig, Proc. Natl. Acad. Sci. 94(1), pp. 148-150 

(1997). 

[15] F. Ritort, C. Bustamante, I. Tinoco, Proc. Natl. Acad. 

Sci. 99(21), pp. 13544-13548 (2002) 
[16] A. M. Berezhkovskii, A. Szabo and J. Weiss, Jour. Chem. 

Phys. 110(18) pp. 9145-9150 (1 999). 
[17] H. Touchette. larXiv:0804.0327V l. 

[18] D. M. Zuckermann and T. B. Woolf, Jour. Chem. Phys 
111(21), pp. 9475-9484 (1999). 



